Random restoration projects

Summary plots

# total random sites to create
tot <- nrow(restdat)

id <- stri_rand_strings(tot, 4)
dts <- range(restdat$date)

# rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
# rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
ext <- bbox(tbpoly)

lon <- runif(10 * tot, ext[1, 1], ext[1, 2])
lat <- runif(10 * tot, ext[2, 1], ext[2, 2])
tmp <- SpatialPoints(cbind(lon, lat), 
                     proj4string = crs(tbpoly)
) %>% 
  .[tbpoly, ] %>% 
  .[sample(1:nrow(.@coords), tot, replace = F), ]

restdat_rnd <- tibble(id) %>% 
  mutate(
    date = sample(seq(dts[1], dts[2]), tot, replace = T),
    top = sample(c('hab', 'wtr'), tot, replace = T)    
  )

reststat_rnd <- tibble(id) %>% 
  mutate(
    lat = tmp$lat, 
    lon = tmp$lon
  )

resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

# map by restoration type
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)

# map by date
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)

# barplot of date counts
toplo <- restall_rnd %>% 
  group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) + 
  geom_bar() + 
  coord_flip() + 
  theme_bw() + 
  theme(
    axis.title.y = element_blank()
  ) +
  scale_y_discrete(expand = c(0, 0))

Distance to restoration sites

wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')

head(wqmtch_rnd)

## # A tibble: 6 x 5
##    stat resgrp   rnk    id     dist
##   <int>  <chr> <int> <chr>    <dbl>
## 1    47    hab     1  06R7 3863.373
## 2    47    hab     2  8o2A 4035.380
## 3    47    hab     3  b5PN 5989.463
## 4    47    hab     4  NRgK 7774.976
## 5    47    hab     5  W146 9705.679
## 6    47    hab     6  gTaE 9821.184

Closest

## 
# plots

# combine lat/lon for the plot
toplo <- wqmtch_rnd %>% 
  left_join(wqstat, by = 'stat') %>% 
  left_join(reststat_rnd, by = 'id') %>% 
  rename(
    `Restoration\ngroup` = resgrp,
    `Distance (dd)` = dist
  )
    
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")

# base map
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
  geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)

# closest
toplo1 <- filter(toplo, rnk %in% 1)

pbase + 
  geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest twenty percent

# closest five percent
fvper <- max(toplo$rnk) %>% 
  `*`(0.2) %>% 
  ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))

pbase + 
  geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest all combinations

# closest all combo
toplo3 <- toplo

pbase + 
  geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Summarizing effects of restoration projects

Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.

salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)

## # A tibble: 6 x 3
## # Groups:   stat [2]
##    stat     cmb     cval
##   <int>   <chr>    <dbl>
## 1     6 hab_aft 23.62066
## 2     6 hab_bef 23.45522
## 3     6 wtr_aft 24.24440
## 4     6 wtr_bef 24.26601
## 5     7 hab_aft 23.80824
## 6     7 hab_bef 24.33082

head(salchg_rnd)

## # A tibble: 6 x 4
##    stat     hab     wtr     cval
##   <int>  <fctr>  <fctr>    <dbl>
## 1     6 hab_aft wtr_aft 23.93253
## 2     6 hab_aft wtr_bef 23.94334
## 3     6 hab_bef wtr_aft 23.84981
## 4     6 hab_bef wtr_bef 23.86062
## 5     7 hab_aft wtr_aft 24.39695
## 6     7 hab_aft wtr_bef 24.11447

Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.

wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)

## # A tibble: 4 x 5
##       hab     wtr              data       crv                    prd
##    <fctr>  <fctr>            <list>    <list>                 <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>

Discretization of salinity conditional probability distributions:

salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd

## # A tibble: 8 x 5
##       hab     wtr      qts       brk  clev
##    <fctr>  <fctr>    <dbl>     <dbl> <dbl>
## 1 hab_aft wtr_aft 25.89943 0.4078706     1
## 2 hab_aft wtr_aft 29.44066 0.8618699     2
## 3 hab_aft wtr_bef 25.59620 0.3693579     1
## 4 hab_aft wtr_bef 29.38185 0.8551371     2
## 5 hab_bef wtr_aft 26.36323 0.4313804     1
## 6 hab_bef wtr_aft 29.94623 0.8684977     2
## 7 hab_bef wtr_bef 26.24658 0.4148461     1
## 8 hab_bef wtr_bef 29.98210 0.8702296     2

A plot showing the breaks:

toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>% 
  unnest
ggplot(toplo, aes(x = cval, y = cumest)) + 
  geom_line() + 
  geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr) +
  theme_bw()

Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.

# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')

# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>% 
  group_by(hab, wtr) %>% 
  nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>% 
  rename(
    salev = cval.y, 
    cval = cval.x
  ) %>% 
  group_by(hab, wtr) %>% 
  nest %>% 
  left_join(salbrk_rnd, by = c('hab', 'wtr')) %>% 
  mutate(
    sallev = pmap(list(data, levs), function(data, levs){
      # browser()
      out <- data %>% 
        mutate(
          saval = salev,
          salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          salev = as.character(salev)
        )
      
      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')

chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')

chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>% 
  print(n = nrow(.))

## # A tibble: 24 x 6
##        hab     wtr salev       qts       brk  clev
##     <fctr>  <fctr> <chr>     <dbl>     <dbl> <dbl>
##  1 hab_aft wtr_aft    lo 12.599474 0.4362990     1
##  2 hab_aft wtr_aft    lo 17.177905 0.8734528     2
##  3 hab_aft wtr_aft    md  6.709566 0.2432915     1
##  4 hab_aft wtr_aft    md  8.343154 0.7019913     2
##  5 hab_aft wtr_aft    hi  4.504318 0.4580072     1
##  6 hab_aft wtr_aft    hi  5.367021 0.8636019     2
##  7 hab_aft wtr_bef    lo 12.123996 0.3361189     1
##  8 hab_aft wtr_bef    lo 16.457887 0.7989985     2
##  9 hab_aft wtr_bef    md  8.037491 0.3559850     1
## 10 hab_aft wtr_bef    md 10.811197 0.8275772     2
## 11 hab_aft wtr_bef    hi  3.968360 0.3206374     1
## 12 hab_aft wtr_bef    hi  4.668884 0.7478114     2
## 13 hab_bef wtr_aft    lo 11.818910 0.3924014     1
## 14 hab_bef wtr_aft    lo 16.028895 0.8215811     2
## 15 hab_bef wtr_aft    md  6.693036 0.4236565     1
## 16 hab_bef wtr_aft    md  8.007187 0.8770172     2
## 17 hab_bef wtr_aft    hi  4.780222 0.4801252     1
## 18 hab_bef wtr_aft    hi  5.800910 0.8786539     2
## 19 hab_bef wtr_bef    lo 11.457732 0.3152265     1
## 20 hab_bef wtr_bef    lo 15.537479 0.7633651     2
## 21 hab_bef wtr_bef    md  8.090232 0.5215649     1
## 22 hab_bef wtr_bef    md 10.561438 0.9094215     2
## 23 hab_bef wtr_bef    hi  4.184045 0.3385945     1
## 24 hab_bef wtr_bef    hi  4.999934 0.7635401     2

Final combinations long format:

chlbar_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  mutate(
    data = map(data, function(x){
      
      brk <- x$brk
      out <- data.frame(
        lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
      )
      
      return(out)
      
    })
  ) %>% 
  unnest %>% 
  gather('chllev', 'chlval', lo:hi) %>% 
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi')),
    chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
  )
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')

chlbar_rnd %>% 
  print(n = nrow(.))

## # A tibble: 36 x 5
##        hab     wtr  salev chllev    chlval
##     <fctr>  <fctr> <fctr> <fctr>     <dbl>
##  1 hab_aft wtr_aft     lo     lo 0.4362990
##  2 hab_aft wtr_aft     md     lo 0.2432915
##  3 hab_aft wtr_aft     hi     lo 0.4580072
##  4 hab_aft wtr_bef     lo     lo 0.3361189
##  5 hab_aft wtr_bef     md     lo 0.3559850
##  6 hab_aft wtr_bef     hi     lo 0.3206374
##  7 hab_bef wtr_aft     lo     lo 0.3924014
##  8 hab_bef wtr_aft     md     lo 0.4236565
##  9 hab_bef wtr_aft     hi     lo 0.4801252
## 10 hab_bef wtr_bef     lo     lo 0.3152265
## 11 hab_bef wtr_bef     md     lo 0.5215649
## 12 hab_bef wtr_bef     hi     lo 0.3385945
## 13 hab_aft wtr_aft     lo     md 0.4371538
## 14 hab_aft wtr_aft     md     md 0.4586999
## 15 hab_aft wtr_aft     hi     md 0.4055948
## 16 hab_aft wtr_bef     lo     md 0.4628796
## 17 hab_aft wtr_bef     md     md 0.4715923
## 18 hab_aft wtr_bef     hi     md 0.4271740
## 19 hab_bef wtr_aft     lo     md 0.4291797
## 20 hab_bef wtr_aft     md     md 0.4533607
## 21 hab_bef wtr_aft     hi     md 0.3985287
## 22 hab_bef wtr_bef     lo     md 0.4481385
## 23 hab_bef wtr_bef     md     md 0.3878566
## 24 hab_bef wtr_bef     hi     md 0.4249455
## 25 hab_aft wtr_aft     lo     hi 0.1265472
## 26 hab_aft wtr_aft     md     hi 0.2980087
## 27 hab_aft wtr_aft     hi     hi 0.1363981
## 28 hab_aft wtr_bef     lo     hi 0.2010015
## 29 hab_aft wtr_bef     md     hi 0.1724228
## 30 hab_aft wtr_bef     hi     hi 0.2521886
## 31 hab_bef wtr_aft     lo     hi 0.1784189
## 32 hab_bef wtr_aft     md     hi 0.1229828
## 33 hab_bef wtr_aft     hi     hi 0.1213461
## 34 hab_bef wtr_bef     lo     hi 0.2366349
## 35 hab_bef wtr_bef     md     hi 0.0905785
## 36 hab_bef wtr_bef     hi     hi 0.2364599

Discretesize chlorophyll data, all stations:

# discretize all chl data by breaks 
chlbrk_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>% 
  mutate(
    lev = pmap(list(data, levs), function(data, levs){
      
      out <- data %>% 
        mutate(
          lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          lev = as.character(lev)
        )

      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest %>% 
  rename(
    chlev = lev, 
    chval = cval
    )

save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')

A bar plot of splits:

ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  facet_grid(hab ~ wtr) +
  theme_bw()

A plot showing the breaks:

toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>% 
  unnest %>%
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi'))
  )
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) + 
  geom_line() + 
  geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr, scales = 'free_x') +
  theme_bw()